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The degradation of this zinc-dialkyl-dithiophosphate (ZDDP), at a temperature less than 246K often leads to the 
release of phosphorus, sulphur, and zinc which are indirectly responsible for the emission of poisonous gas from the 
exhaust pipe of the motor cars. Four QSPR mathematical models were generated from 39 structures of lubricant 
additives (LAs) and the structural features were found to corresponds to the coefficients of; internal correlation (Rô) 
of 0.95, adjusted squared correlation (Radj) of 0.94, Cross-validation (Q2cv) of 0.90, and the external validation 
(R*pred) of 0.54. The model suggests that new LAs with improved onset temperatures (Tonset) Could be designed by 
interpreting and increasing the value of the molecular descriptor such as IC5 (Information Content 
index/neighborhood symmetry of 5-order) and Ve (V total size index/weighted by Sanderson electronegativity) and at 
the same time decreasing the values of RDF080m (Radial Distribution Function-080/weighted by mass), RDF110m 
(Radial Distribution Function-110/weighted by mass), P2v (2nd component shape directional WHIM index/weighted 
by Van der Waals volume) and R1e+ (R maximal autocorrelation of lag 1/weighted by Sanderson electronegativity). 
Moreover, the LAs with an experimental onset temperature of 351.6K agreed with the predicted onset temperature of 
351.7K13a. And was also in agreement with the result of molecular dynamics simulations in which the LAs with the 
best dynamic binding energy of -2112.06 kcal/mol was tightly bounded on the simulated DLC mechanical coated 
boundary inter-surface and was also found to be better than the commercial LAs, ZDDP in term of binding energy 
and onset temperature. This investigation will help in rational additive design and synthesis of new and better 
selective Las 


Keywords: Lubricant Additives (Las); Diamond-Like Carbon (DLC); DFT; GFA; Molecular Dynamic (MD); 
simulations; Quantitative Structure-Properties Relationships (QSPR) 


Introduction 


Research showed that energy loss from ineffective boundary lubrication of mechanical components is accounting for 25% of 
important energy needed while 15% of the loss in internal energy is due to mechanical frictional loss (Yanli, (2017), and Becker, 
(2004)). Innovative machine technology has made it clear that the requests put upon each viable oil were not met by even the 
monetarily accessible lubricants. To avoid friction and machine wearing, researchers are searching for reliable and less costly means 
to designs different types of lubricant additives (LAs) that could increase the known properties of lubricants and will not decompose 
easily when blended with the base stock at higher machines lubricating dynamic temperature (Kato, (2011); El-Eisawy et al., (2015); 
Ahmad et al., (2018); Podgornik, (2001); Moriguchi et al., (2014)). LAs are chemical organic compounds used to increase the 
performance of base oil stock (Ahmad et al., (2018); Podgornik, (2001)). LAs help to prolong the machine engine lifespan by 
providing the necessary performance for smooth and better operation (Kato, (2011); El-Eisawy et al., (2015); Ahmad et al., (2018); 
Podgornik, (2001)). Scientists have been effectively performing different researches on the association of different DLC boundary 
coatings with various kinds of LAs to acquire better machine tribology understanding of boundary interfaces (Neville et al., (2007); 
Kano et al., (2005); (2005), Kalin and Velkavrh, (2013); Aboulthana et al., (2019); Forsberg et al., (2013); Chiro et al., (2010)). LAs 
such as Pyrimidine, natural acids, and pyrazines have emerged following the used zinc-dialkyl-dithiophosphate (ZDDP) as 
multifunctional LAs. This widely used ZDDP was reported to hurt the environment and was also reported to decompose at a 
temperature of less than 246K (Chiro et al., (2010); Zhongyi et al., (2020)). Therefore, urgent need to look for better and ecologically 
friendly LAs that can withstand diverse high lubricating working temperatures cannot be overemphasized. Quantitative Structure- 
Properties Relationships (QSPR) and molecular dynamics simulations are computational methods widely used as a research tool 
(Kubinyi, (1997); Robinson et al., (1999); Ivanciuc et al., (2000); Wong et al., (2002)). QSPR relates molecular properties with its 
molecular structure of interest. The advantage of the QSPR technique is the capacity to figure the properties of new hypothetical 
chemical substance and it serves resources and gives direction before the real experimental study (Abdulfatai et al., (2019a); 
Mohammad and Zohreh, (2013); Abdulfatai et al., (2019b)). 
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Moreover, the molecular dynamics simulations main success was the possibility to calculate the dynamic binding strength (binding 
energy) of a molecule on interfaces accurately before carrying out the complicated expensive experimental study (Abdulfatai et al., 
(2019a); Mohammad and Zohreh, (2013); Abdulfatai et al., (2019b)). This investigation was aimed at designing better LAs that can 
resist high dynamic boundary temperatures before it decomposes and to determine the dynamic binding strength of LAs in the 


boundary sliding interfaces via QSPR and MD. 


1 Materials and Methods 
1.1 Data sets, Molecular Descriptors Generation, and model building 


The 2D structure of LAs (Table 1) with their onset temperatures were retrieved from the literature and used for this research (Deng, 


(2009); Jeng, (2009); Jiao, (2011)). 


Table 1 Experimental, Normalized and Predictive Lubricant Additives Onset Temperature (Tonset ) 
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The 2D structures of all the 39 LAs were drawn with ChemDraw software and later transferred and converted to 3D by Spartan’ 14 
version 1.1.2 software. Using molecular mechanics (MMFF), these 3D structures of all the LAs were then subjected to energy 
minimization (Hashem, (2011)) after which complete geometry optimization was carried out with Density Functional Theory 
(B3LYP/6-311++ G**). The optimized LAs were transferred to Dragon 6.0 software toolkits (Todeschini et al., (2017)) where about 
3556 molecular descriptors were generated. The onset temperature of these LA compounds was measured as T4,4,(K) were 
normalized and further expressed as p'""*', Using a random selection method, all the LAs were grouped into 25 and 14 sets. Material 
studio version 8.0 software was utilized to carry out correlation analysis and development of the QSPR models with the genetic 
function algorithm method (GFA). GFA as a method, as the ability to incorporate the idea of natural evolution (Rogers and 
Hopfinger, (1999); Ching-Wen et al., (2014)). The quality assurance of the generated QSPR mathematical models was accessed by 
internal (R, Q Ra) and predictive i prea-) validation parameters (Equations 1-4) in comparison with the acceptable QSPR 
parameters (Ravinchandran et al., (2011)). The relative relevance of every generated descriptor contribution to the developed model 
was calculated in-term of mean effect (MF) value in equation 5 (Ravinchandran et al., (2011)). 
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Where Y y, -observed/experimental training set properties, 
Y ic-calculated/predicted training set properties, Y obstet) 
observed/experimental properties for the test set, Y 44,4, = predicted test set 
properties, N= number of molecules in the data, Yp =represent the predicted 
properties of the training set, Y=observed properties of the training set, Y m 
the mean activity value of the training set, R?-determination coefficient, 
p=number of descriptors in the QSPR model, N-1-p- degree of freedom, j= 
model’s generated descriptor, Bj = descriptor's coefficient, dij= training 
set's data matrix descriptor value, m = sum of model’s descriptors while n = 
summation of training set's molecules. 


Fig. 1 3D structures of DLC substrate and LAs 
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1.2 Molecular dynamics (MD) simulation research 


To understand, access the structural and dynamic binding energy between the DLC simulated coated interfacial sliding surfaces and 
potentially designed lubricant anti-oxidant additives, molecular dynamics (MD) simulation investigation was studied (Wymyslowski 
et al., (2008)). The 3D structure (Figure 1) of hydrogen-containing DLC crystal surface which was reported in many works of 
literature (Wymyslowski et al., (2008)) to be better allotropy of carbon in term of wearing resistance ability and antioxidant in the 
sliding interface was constructed from the carbon (C) model using Materials studio 8.0 simulation software. This was done by 
cleaving the carbon surface at point 1.1 .0 (Height, Breadth Length), Top (1.006 A), and Thickness (24.121 A) into the crystal unit. 
Moreover, the repeated units (supercell) of the carbon crystal units were formed at (Height) U (9) and (Breadth) and V (8) while 
hydrogenation of the supercell, vacuum slab, and geometric optimization was performed. Similar steps were also followed to build 
steel. The 2D structures of the lubricant antioxidant additives were drawn with Chemdraw software and were then converted to 3D 
structures by materials studio version 8.0 simulations Program, then optimized and saved by the same software (Fig. 1). COMPASS 
II which is a strong and well-created power field (than COMPASS) that was determined dependent on the fitting capacity against a 
wide scope of inorganic and synthetic organic compounds (Wymyslowski et al., 2008) was selected in the Materials studio 8.0 
Software. Molecular Dynamic (MD) simulation binding energy (B.E) calculations were done subsequent after presenting the 
minimized additive compound into the simulation vacuum slab of optimized steel and DLC crystal (24.82 A x24.82 A x45.27 A) 
surface at 350.15 K and over a range of inter-surface separations. The interaction binding energy of all the designed additives on the 
diamond like carbon (DLC), was determined by using equation 6 (Zhao et al., (2014)). 


E Dynamic Binding Energy =E tora (E LubricantAdaitivet Epic Surface) (6) 


2 Results and discussion 
2.1 Analysis of QSPR 


After using the Genetic Function Algorithm (GFA) method, GFA in material studio software to carefully developed the QSPR model 
from 39 structures of LAs and the onset properties, four QSPR mathematical models were generated. All these generated 
mathematical models along with their statistical correlation coefficients were subsequently subjected to QSPR validation parameters 
(Ravinchandran et al., (2011)). The first, second, third and fourth models were found to have coefficients of; internal correlation (Rô) 
of 0.95, 0.94, 0.94, 0.94, adjusted squared correlation (Radj) of 0.40, 0.93, 0.93, 0.92, Cross-validation (Q)) of 0.90, 0.90, 0.80, 0.89, 
and the external validation (R’pred) of 0.54, 0.52, 0.52, 0.51, respectively. But after careful comparison, the 1st model was revealed 
to better comply with the standard validation requirements (Ravinchandran et al., 2011). 


Model 1: 


ponse) = 0.07*1C5 — 0.02*RDF080m — 0.09*RDF110m — 1.23*p2v + 0.00002*Ve — 3.22*R1e + +2.56.R2 = 
0.95, RZext = 0.54, R2adj = 0.94,Q2cv = 0.90. 


Model 2: 
pense) = Q.08'IC5 — 0.02*RDF080m — 0.009*RDF0800m — 1.27*RDF110 + 0.00002*P2v — 3.28*Rle + +2.55.R2 
= 0.94,R2ext = 0.52, R2adj = 0.93,Q2cv = 0.90. 
Model 3: 
pense) = 0,12*IC5 — 0.03*RDF080m — 0.03*Mor1le + 0.38*Dv + 0.18' H3p — 2.47*Rle + +1.91. R2 = 0.94,R2ext 
= 0.52, R2adj = 0.93,Q2cv = 0.90 
Model 4: 


ponse) = 0.08*IC5 — 0.024*RDF080m — 0.012*RDF110 — 1.32*P2v + 0.00002*Vm — 3.325*R1e + +2.56.R2 
= 0.94,R2ext = 0.51,R2adj = 0.92,Q2cv = 0.89. 
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The confirmation assessment further proved that the 1* model’s internal coefficient prediction (R^) of 0.95 value conformed with R? 
coefficient (training set graph) of 0.95 (Figure 2) and the external coefficient of the predicted (ge of 0.54 value was also in 
conformity with R? (test set graph) value of 0.54 in the plotted 
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Fig. 2 Plot of training set 
Fig. 3 Plot of test sets 


graph (Figure 3). Moreover, there was no problem of multi-colinearity since model 1’s descriptors coefficient correlation was very 
low (Table 2) (David et al., (2004); Costa et al., (2014)). Therefore, in QSPR model 1, the positive coefficient signs in the molecular 
descriptors such as IC5 and Ve while the negative coefficient sign of RDF080m, RDF110m, P2v, and RJe+ implies that increase in 
former and decrease in latter descriptors will increase the onset temperatures (7 nse) properties at which the lubricant additives 
decompose during the mechanical boundary lubrication. Since all the required statistical analysis carried out have proved that most of 
the 39 LAs used in this study can withstand high temperature without being decomposing easily during mechanical boundary 
lubrication of machine of interest. It is also important to know those that can enhance the properties of base oil at high onset 
temperature. Applicability domain (AD) was a help in this case. AD statistical analysis was performed to determine the set of 
molecular additives found within the domain and the outliers/influential additives compounds outside the domain or danger 
boundary, h*(Eqn. 7) (Netzeva et al., (2005)). Where m=number of descriptors that appear in a QSPR linear model, n=number of 
training set molecule only. 


. 3(m+1) Table 2 Model I descriptors inter-correlation matrix 


h» n V) Descriptor ICS RDF080m RDF110m P2v Ve Rle+ 

IC5 1 

A ee look at Figure 4 shows that all the ree wee RDFO80m 023 1 

found within the domain of the graph except additive with 2 

serial number 39 (Table 1) that was noticed to be outside RDF110m 0.19 0.45 1 

the danger limit (h*) of 0.83. Therefore, this compound is — jj, 0.59 0.48 -0.007 1 

called outlier and it cannot be used as a template to design - 

other LAs because it is not reliable and its property may Ve 0.30 0.42 0.93 0.12 1 

be gotten in error during the experimental characterization pjey 028 -0.73 -0.56 0.50 0.52 1 


(Netzeva et al., (2005)). 
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2.2 Molecular Dynamic Simulation Analysis 


To determine the dynamic binding energies of all the LAs under investigation, computational MD simulation studies were carried out 
between LAs and the boundary DLC coated surface. The binding energy scores in Table 3 show the individual binding energies of 
lubricant additives, DLC, total and binding energies of the complex. All the lubricating oil additive compounds (Table 3) were found 
to bind strongly with the DLC mechanical sliding interface. Besides, LAs with serial number 13a and with a dynamic binding energy 
of -2234.05 kcal/mol was found to show excellent binding energy than other co-lubricant additives. This study revealed that this LA 
with serial number 13a could form a stable film on the surface of the DLC boundary coated surface than other co-additives. Figure 5 
shows the best LA- DLC complex with the best binding energy. This binding energy of -2234.05kcal/mol generated was found to be 
better than the one reported by Costa and his co-workers in 2011 (Costa et al., (2011)) and it was also better than the dynamic binding 
energy of ZDDP (Table 3). 


Table 3 Lubricant Additives Binding Energies 


Comp. No Ei Additives (kcal/mol) Epic (kcal/mol) Exot (kcal/mol) EBinding Energy (kcal/mol) 
la -35.29 -970.04 -1024.3 -1959.06 
2a 41.46 -970.04 -940.75 -1952.25 
3a 56.30 -970.04 -933.49 -1959.83 
4b 56.42 -970.04 -937.16 -1963.62 
5b 91.98 -970.04 -893.3 -1955.32 
6b 91.49 -970.04 -788.57 -1850.09 
7b 91.59 -970.04 1710.13 648.50 
8b 91.20 -970.04 -832.65 -1893.89 
9a 238.43 -970.04 -785 -1993.47 
10b -4.81 -970.04 -94.29 -1059.51 
11b -4.62 -970.04 36.47 -928.947 
12 -4.63 -970.04 -74.51 -1039.91 
13a 240.21 -970.04 -1023.8 -2234.05 
14a -26.00 -970.04 605.39 -338.65 
15b 36.90 -970.04 -951.42 -1958.36 
16a 33.5 -970.04 832.03 -171.50 
17b 23.95 -970.04 3443.24 2449.25 
18b 29.23 -970.04 3285.99 2286.71 
19b 35.46 -970.04 722.96 -282.54 
20b 16.10 -970.04 -972.93 -1959.07 
2la 14.31 -970.04 -1011.7 -1996.09 
22b 13.95 -970.04 -965.21 -1949.2 
23b 54.28 -970.04 -929.46 -1953.77 
24b 56.79 -970.04 -915.82 -1942.64 
25a 52.14 -970.04 -923.42 -1945.61 
26b 53.46 -970.04 -921.24 -1944.73 
27a 83.19 -970.04 -652.05 -1705.29 
28b 144.34 -970.04 853.375 -261 

29b 195.17 -970.04 -742.07 -1907.28 
30b 104.80 -970.04 7487.52 6412.68 
31b 94.97 -970.04 3557.07 2492.06 
32b 2.99 -970.04 6379.47 5406.43 
33a 275.67 -970.04 2736.1 1490.40 
34b -9.43 -970.04 7169.77 6809.17 
35b 34.90 -970.04 -969.17 -1974.12 
36a 20.14 -970.04 -973.15 -1963.33 
37a 33.22 -970.04 -983.84 - 1987.09 
38b 31.17 -970.04 1640.3 639.09 
39a 36.80 -970.04 727.42 -279.41 
ZDDP 262.62 -970.04 18502.31 24365.40 
EL. aaditives=Lubricant Additives Epic = Diamond Like-Carbon, Erota = Total Energy of the complex Training set: “Test set 
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Conclusions 


QSPR's principle which was based on chemical properties-structure relationship and MD simulation which was based on the 
calculation of dynamic binding energy of compound of interest were used to carry out the onset temperature of LAs study. Four 
QSPR mathematical models were generated from 39 structures of LAs and the onset properties. The 1* model was found to better 
comply with the standard validation requirements. The generated physicochemical descriptors from the model correspond to the 
structural features and were found to have coefficients of; internal correlation (R°) of 0.95, adjusted squared correlation (R? adj) of 
0.94, Cross-validation (Qa) of 0.90, and the external validation (R’pred) of 0.54. The model suggests that new LAs with improved 
onset temperatures (T nse.) could be designed by interpreting and increasing the value of the molecular descriptor such as IC5 and Ve 
and at the same time decreasing the values of RDF080m, RDF110m, P2v, and R/e+. Moreover, the LAs with an experimental onset 
temperature of 3 51.6 K agreed with the predicted onset temperature of 351.7 K (13a, Table 1). And was also in agreement with the 
result of molecular dynamics simulations in which lubricating oil additive number 13a, Table 3, with the best dynamic binding 
energy of -2112.06 kcal/mol. These two methods could be adopted by researchers and industrialists/engineers to design an improved 
lubricant additive property before wet experiments since they have proven to be good research tools due to the structure-properties 
correlation capability. 


Nomenclature 

QSPR =Structure-Properties Relationships [-] 

DLC =Diamond like carbon [-] 

DL Cubricant additive =Energy of the diamond like carbon [-] 

Dv =D total accessibility index / weighted by Van der Waals volume [-] 

E =Energy [Kcal/mole] 
Lubricant additive =Energy of the lubricant additive [Kcal/mole] 
Exotal =Total Energy [Kcal/mole] 
GFA =Genetic function algorithm [-] 

ICS =Information Content index/neighborhood symmetry of 5-order [-] 

LAs =Lubricant additives [-] 

MD =Molecular dynamic [-] 

P2v =2nd component shape directional WHIM index / weighted by Van der Waals volume [-] 

Q2cv =Cross-validation [-] 

RDF =Radial Distribution Function [A?] 
RDF080m =Radial Distribution Function - 080 / weighted by mass [-] 
RDF110m =Radial Distribution Function - 110 / weighted by mass [-] 

Rle+ =(R maximal autocorrelation of lag 1 / weighted by Sanderson electro negativity [-] 

Rle =R autocorrelation of lag 1 / weighted by Sanderson electro negativity [-] 

Tonset =Onset temperature [-] 

Ve =V total size index / weighted by Sanderson electro negativity [-] 

Vm =V total size index / weighted by mass [-] 
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